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APPLICATION OF MONTE CARLO TECHNIQUES TO OPTIMIZATION 
OF HIGH-ENERGY BEAM TRANSPORT IN 
A STOCHASTIC ENVIRONMENT 

By Russell V. Parrish, James E. Dieudonne, 
and Tassos A. Filippas 
Langley Research Center 

SUMMARY 

An algorithm employing a modified sequential random perturbation, or creeping 
random search, has been applied to the problem of optimizing the parameters of a high- 
energy beam transport system. The stochastic solution of the mathematical model for 
first-order magnetic-field expansion allows the inclusion of state -variable constraints, 
and the inclusion of parameter constraints allowed by the method of algorithm application 
eliminates the possibility of infeasible solutions. The mathematical model and the algo- 
rithm were programed for a real-time simulation facility; thus, two important features 
are provided to the beam designer: (1) a strong degree of man-machine communication 

(even to the extent of bypassing the algorithm and applying analog- matching techniques), 
and (2) extensive graphics for displaying information concerning both algorithm operation 
and transport -system behavior. Chromatic aberration was also included in the mathe- 
matical model and in the optimization process. 

Results presented show this method as yielding better solutions (in terms of resolu- 
tions) to the particular problem than those of a standard analog program and well as dem- 
onstrating flexibility, in terms of elements, constraints, and chromatic aberration, allowed 
by user interaction with both the algorithm and the stochastic model. Also presented are 
examples of slit usage and a limited comparison of predicted results and actual results 
obtained with the 600 MeV cyclotron at the NASA Space Radiation Effects Laboratory. 

INTRODUCTION 

With the advent of large-accelerator projects over the past 15 years, a demand has 
arisen for optimization of not only the design of the accelerator but also the design of the 
beam-handling systems that must necessarily be coupled with the accelerator. As pointed 
out in reference 1, the high cost of producing the strong magnetic fields over large vol- 
umes necessary for high-energy beam optics makes optimization a necessity. 



Past efforts in the area of beam -transport optimization have generally fallen into 
two classes of techniques: ’’analog matching" with a high-speed repetitive analog com- 
puter and iterative deterministic methods for parameter optimization implemented on a 
general-purpose digital computer. References 1 to 3 describe the general methods used 
in the technique of analog matching: Generally in designing transport systems the results 
of parameter variations in the system are continuously displayed on an oscilloscope as 
the designer seeks by trial and error to determine the focusing strengths and separation 
distances necessary to give the desired results. Aside from the problems of analog res- 
olution and repeatability, which are usually present even with advanced equipment, the 
major disadvantage with this technique is that it relies heavily on the experience of the 
analyst, especially as the number of parameters increases. It is extremely time consum- 
ing and is usually based on visual determination of improvement; thus, the problem of 
visual resolution arises. 

However, as references 1 to 3 strongly emphasize, the man-machine communication 
available with the analog computer is often felt to be the overriding consideration when 
choosing between the two classes of techniques available to solve the problem of beam- 
transport design. This man-machine communication is generally sacrificed when the 
analyst chooses to use a digital-computer optimization program. This sacrifice is not 
unwarranted, however, according to references 4 and 5, which emphasize the advantages of 
optimization algorithms over analog methods. These advantages include a performance 
measure for evaluating a given set of parameters, an algorithm that automatically seeks 
an extremum solution of the performance measure (which may or may not be feasible), 
and digital resolution. 

Thus a survey of the literature available in the field reveals strong advocates for 
each class of techniques. The present paper will describe a method that evolved in the 
process of solving a particular problem in beam -transport design. This method combines 
some of the best features of both classes of techniques, that is, a digital algorithm with its 
ensuing resolution, and man-machine communication (complete communication even to the 
extent of bypassing the algorithm and applying analog- matching techniques). In addition, 
the method allows for the inclusion of state -variable constraints and parameter con- 
straints in the mathematical model. Other features of the program include stochastic 
representation of initial conditions and chromatic aberration and the introduction of the 
aberration during the optimization process. These features make the simulation of the 
beam transport system more realistic than in most previous work (e.g., ref. 1). 

The method will be presented in the context of solving a specific problem in beam- 
transport design, namely, that of transporting a beam of charged particles from one given 
area to another through the use of four quadrupole magnets (other configurations of ele- 
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ments are allowable). In tins problem five parameters are available for the optimization 
process. The maximum number of parameters that can be varied at one time is deter- 
mined by the choice of the optimization algorithm and the amount of computing time nec 
essary to evaluate the performance measure. However, the method may be used to opti- 
mize several elements at a time along the transport system, with the resulting beam used 
as input to the next set of elements in large systems. 

SYMBOLS 

ajjbj constraint constants, meter"! 

C value of cost function with present parameter values, dimensionless 

Cp value of cost function at last success, dimensionless 

C x cost function corresponding to horizontal axis, dimensionless 

C z cost function corresponding to vertical axis, dimensionless 

D2 separation distance between second and third quadrupoles, meters 

e electron charge, coulombs 

F ca i calculated test statistic, dimensionless 

F gg^n,np'j F value for 95-percent confidence level with ^n,n^j degrees of freedom, 
dimensionless 

ith magnetic-field gradient, webers/meter^ 
alternate hypothesis for a statistical test, dimensionless 
null hypothesis for a statistical test, dimensionless 
magnetic focusing strength of ith quadrupole, meter"! 
effective field width, meters 

number of rays which passed through the system, dimensionless 


Si 

h a 

H 0 

Ki 

h 
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n 

n P 

N 

P f,X 

P f,z 

min 

( P f> x )min 

( P f> x )max 

j z ) min 

(P f z ) 

V 1 > z ^max 

Po 

Ap 

r 

s 

As 

s e,x 

s e,z 

s f 



number of rays transmitted with present parameter values, 
dimensionless 

number of rays transmitted at last success, dimensionless 

percentage of rays which did not pass through the system, 

n ~ m x 100, dimensionless 
n ’ 

horizontal focal point, meters 
vertical focal point, meters 

minimum value allowed for distance between horizontal and vertical focal 
points of the first doublet, meters 

minimum value allowed for horizontal focal point of the first doublet, meters 
maximum value allowed for horizontal focal point of the first doublet, meters 
minimum value allowed for vertical focal point of the first doublet, meters 
maximum value allowed for vertical focal point of the first doublet, meters 
average momentum, kilogram -meters/second 
momentum perturbation, kilogram -meters/second 
thickness of slit, meters 

displacement along focal axis of transport system, meters 
s out - s n> meters 

position along focal axis of extremum in x, meters 
position along focal axis of extremum in z, meters 
position along focal axis at which cost function is evaluated, meters 
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s h,B 

s h,E 

s n 

s out 

S V,B 

s v,E 

*cal 

to 

*1 

t.05( n + n p 

Wdif 

w N 

x(s) 

*d 

F( Sf ) 

x sl 

z(s) 


beginning position of horizontal slit, meters 
ending position of horizontal slit, meters 

position of entrance to field (magnetic or drift) along focal axis, meters 

position of exit from field (magnetic or drift) along focal axis, meters 

beginning position of vertical slit, meters 

ending position of vertical slit, meters 

calculated test statistic, dimensionless 

t value associated with type I error, dimensionless 

t value associated with type n error, dimensionless 

- 2^ Student's "t" for 95-percent confidence level with (n + iip - 2j degrees 
of freedom, dimensionless 

weight corresponding to difference between cost function for vertical plane 
and cost function for horizontal plane, dimensionless 

weight corresponding to percentage of beam which passed through the system, 
dimensionless 

displacement along horizontal axis of transport system, meters 

desired value of horizontal displacement at end point, meters 

mean value of the horizontal slope at the end point of all rays of the beam 
which passed through the system, dimensionless 

displacement of beam in horizontal plane parallel to focal axis, used to deter- 
mine slit opening, meters 

displacement along vertical axis of transport system, meters 
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a ; 


a 2 




cr 


co. 


X 


COn 


Subscript: 


desired value of vertical displacement at end point, meters 

mean value of the vertical slope at the end point of all rays of the beam 
which passed through the system, dimensionless 

aperture radius, meters 

detectable difference desired, dimensionless 

variance of cost function for present parameter values, dimensionless 
estimate of cost-function variance at last success, dimensionless 

variance of cost function at last success, dimensionless 

estimate of cost-function variance for present parameter values, 
dimensionless 

slit width along x-axis, meters 
slit width along z-axis, meters 


d desired value 

Primes indicate derivatives with respect to the focal axis (s-axis). 


MATHEMATICAL MODEL 


The basic beam-transport problem discussed in this paper is the design of a high- 
energy beam transport system consisting of four quadrupoles, with variable magnetic 
focusing strengths and separation distances, and slits with variable width, thickness, and 
separation distance. This system will transform a charged-particle beam of known initial 
conditions into a beam meeting some desired phase -plane configuration defined by a math- 
ematical formula known as the cost function. In the design of this particular system, 
illustrated in figure 1 , five system parameters were to be determined. These were the 
four magnetic focusing strengths /Kj, K2 ? K3, and K4) and the separation distance D2 
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between the doublets (pairs of quadrupoles). All other system parameters are considered 
fixed. In general, the method used to determine the parameters imposes no restrictions 
as to number or types of the parameters. 

For applications of interest here, the differential equations of motion for traversing 
a quadrupole field (ref. 1, p. 9) 



can be reduced to the uncoupled, linear equations (refs. 1 and 6) 


x" + K^x = 0 


z" - iqz = 0 


( 2 ) 


where iq, a constant (i = 1,2, 3, 4), is the focusing strength of the field. It represents a 
rectangular, first-order, field model Kq = eg^p 0 (where the ith magnetic potential is 
equal to -g^ • xz) in which focusing occurs in the xs-plane and defocusing in the zs-plane. 
The algebraic solutions of equations (2) are shown in the following equations: 


x ( s out) =x(s n )cosNK i As] + 


x'(s n ]sin (fiq as) 


li'Ki 


(3 a) 


x ’( s out) = - x ( s n) jiq sin(jiq As) + x'fs n )cos (jiq As) 


z^s n ] sinh ^\(Ki 


As 


Ki 


z ( s out) = z ( s n)cosh^|iq As) 


z'^Sout) = z ^ s n)\j K i sinh^jiq As) + z'^s n )cosh^]iq As 


(3b) 

(3c) 

(3d) 


As - s out - s n 
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The inverse case, namely, focusing in the zs-plane and defocusing in the xs-plane, may 
be implemented by using equations (2) as shown, but with the sign of Kj changed. 

The drift-space translational equations 



also have algebraic solutions: 

x ( s out) = x ( s n) + x '( s n) As (5a) 

x '( s out) = x ’( s n) (5b) 

z(sout) = z(s n ) + z'^Sj^As (5c) 

z '( s out) = z '(s n ) (5d) 


Throughout the entire translation along the s-axis, the position variables are sub- 
ject to the physical constraints imposed by the aperture size of each successive 
quadrupole: 


x2(s) + z^(s) 5 a? 


(6) 


where a ^ are the aperture radii for the four quadrupoles (i = 1,2, 3, 4). This constraint 
was imposed at the entrance and exit of each quadrupole as well as within the quadrupole. 
The locations of the extremum values of x and z ^s e x and s e within the focusing 

quadrupole are determined by 


"\ 



(7) 
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where 



r 

0 

( S e,* « °) 


f 

0 

( S e, z •= °) 

s e,x - < 

s* 

epc 

(° S s c',x s As ) 

s e,z - j 

s* 

e ?z 

(° 5 s e,z s As ) 


As 

{ S U > is ) 


As 

( s e,z > As ) 


k. V. 


Again, the equations presented are for focusing in the xs-plane, and a change in sign of 
Kj is necessary for focusing in the zs-plane. 

In addition to the position constraints of equation (6), two types of parameter con- 
straints may be present in the problem. The first type 


-bi =§ Ki =§ bi 

a l = D 2 = a 2 


(i = 1,2, 3 ,4) 

> ( 8 ) 


J 


is imposed by the physical properties of the quadrupoles to be used in the actual experi- 
ment, limiting the values of Kj to feasible ranges and the separation distance between 
the doublets to values meeting the requirements of the experimental laboratory 
dimensions. 

The second type of parameter constraint is introduced by the desire to use slits to 
improve the beam quality. Slits are usually to be introduced at some focal point in either 
the xs- and zs -planes. It is desirable that the focal points fall within the drift space 
between the elements and still allow enough room for slit placement. It is also desirable 
that the focal points be separated from one another by sufficient distance for placement 
of slits of different thickness at each focal point. These constraints are represented by 


M min * p £ ,x s Kx) 

( Pf ’ z )min S S M 


'f,x " P f,z 


2 ( ip f) 


min 


max 


max 


(9a) 

(9b) 

(9c) 
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and the focal points are determined by solving equations (3) with initial conditions 
x' = z' =0, and x = z = Constant through the desired elements. The focal points are 
then determined by the points of intersection with the s-axis. 

Slits are implemented as position constraints: 

l x H < X ( s h,B S s S s 1,,e) 

l z < s >l < X ( S V,B S S S S v,e) 

assuming no penetration (complete absorption), and can be introduced separately, jointly, 
or not at all. With the above constraints imposed on the system, the problem is to find 
a set of parameters which when placed in the system will yield a "best fit" to some pre- 
determined phase-space curve defined in the cost function. The algorithm chosen to per- 
form this optimization is explained in the section of this paper entitled "Monte Carlo 
Method. " 

One primary aberration, chromatic aberration, was taken into account in the mathe- 
matical model. This aberration arises from the different values of momentum possessed 
by the individual rays. In the equations presented previously, the momenta of the parti- 
cles are considered constant and equal throughout the transport system, that is, 

Kj = egj/p 0 . Usually an optimal case is obtained without the presence of chromatic aber- 
ration, and then the effect of the aberration is examined in the following manner: The 

eg t 

varying momentum is modeled as Kj = , where Ap is a positive or negative 

P 0 +AP 

momentum perturbation. In the present paper, each ray of the beam is assigned a momen- 
tum perturbation based on Gaussian noise. The perturbation is then applied to each focus- 
ing strength. In this manner, the momentum of an individual particle of the beam is con- 
stant throughout the focusing system, while individual momenta within the beam vary. 
Optimization then takes place in the presence of chromatic aberration. 

ANALOG SIMULATION 

During the preliminary investigation of the problem, the equations of motion 
(eqs. (2)) and the constraint equations (eqs. (8)), were programed on both the GPS 10000 
repetitive analog computer and the EAI 680 analog-hybrid computer. The function Kj, 
which represents the system parameters, was generated on both machines by the use of 
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electronic comparators and available logic components. The variation of system param- 
eters was achieved by tandem potentiometers. Repetitive solutions of the differential 
equations at a rate of 10 runs per second allowed oscilloscope displays of the trajectories 
z(s) and x(s) (beam traces) and phase -plane plots in both planes at any point in the 
system. 

A visual-resolution problem was encountered in the use of this method in that the 
operator found many promising solutions, or sets of system parameters, from which by 
visual means a "best case" could not be determined. To solve this problem, a digital 
program was written for the solution of the differential equations with the predetermined 
sets of system parameters. The digital program, which also served as an independent 
check, solved the resolution problem, but uncovered the fact that the analog-determined 
system parameters did not fulfill the desired phase-plane end conditions. This fact will 
be demonstrated in the "Results" section of this paper. 

REAL-TIME SIMULATION SYSTEM 

The inability of the analog program in meeting required phase -plane end conditions 
led to the programing of the problem on the CDC 6600 Real Time Simulation (RTS) 
Subsystem at the Langley Research Center. This system, which is described in refer- 
ence 7, provides both sufficient accuracy for the problem and man-machine communica- 
tion. In addition to the main computers, the RTS system contains program control sta- 
tions, an analog-to-digital and discrete input system, a digital-to-analog and discrete 
output system, a CDC 250 Series CRT display system, and peripheral units (time-history 
recorders, x-y plotters, card readers, etc). The RTS software provides many "analog- 
like" features such as mode controls, subroutines for displaying and changing parameters, 
and function-sense switches. 

The application of the RTS system to the beam -transport problem allows the use of 
an optimization algorithm to determine an optimal set of system parameters while the 
operator maintains man-machine communication. 

This communication link to the machine allows quick, effective use of the designer's 
knowledge of the system. It enables him to examine different areas of the parameter 
space and determine the merits of particular parameter sets as to likely starting points 
for the optimization algorithm, thereby reducing the time required for convergence. The 
designer also has many options available through subroutines, function-sense switches, 
and the RTS software from which he may change the model of the transport system, 
position-variable and parameter constraints, the form of the desired output of the system, 
and the form of the probability distribution that dictates the initial conditions of the beam 
prior to magnetic-field entrance. 
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Monte Carlo Method 


The individual rays of the incident beam in the model are chosen to represent the 
physical beam for which the transport system is being designed. This is done by choos- 
ing a set of random initial conditions that satisfies the physical requirements. Equa- 
tions (3) are then evaluated to provide the values of the position variables at the entrance 
and exit of each magnetic field, and also the extremum values within the magnetic fields 
are evaluated (eqs. (7)). Any violations of the position-variable constraints are recorded 
(solution of the equations of motion for an individual ray is terminated upon the occur- 
rence of a violation) and may be included in a weighted performance measure with the 
deviation of the transmitted particles from the desired end conditions. This process is 
then repeated and the deviations are accumulated to form the mean performance measure, 
or mean cost function, over the desired sample size. The parameters of the system (four 
focusing strengths and one separation distance) are then perturbed by the optimization 
algorithm and the entire sequence is iterated until convergence is achieved. 

Monte Carlo Algorithm 

The optimization algorithm used in conjunction with the mathematical model was a 
variation of the sequential random perturbation, or creeping random search, suggested by 
Mitchell in reference 8. This algorithm was selected over deterministic methods that 
are gradient dependent, such as steepest descent, conjugate gradient, Newton-Raphson, 
and quasilinearization (see ref. 9) , because of the presence of constraints and a stochastic 
model. The algorithm is described in detail in reference 10. 

A stochastic representation of the initial beam (beam prior to entry into the trans- 
port system) was advantageous in making the simulation of the system more realistic. 

By considering the beam as a conglomerate of rays, each having individual characteristics 
of displacement, slope, and momentum, position-variable constraints may be implemented, 
and the percentage of the beam striking the elements of the transport system may be cal- 
culated. The stochastic model also allows for realistic representation of chromatic aber- 
ration since each ray of the beam instead of the beam as a whole is assigned a momentum 
perturbation. 

The basic algorithm, devoid of any strategies, is best illustrated with the two- 
parameter contour plot shown in figure 2(a). In this simple problem, constant values of 
the cost function are represented as concentric circles in the two-parameter plane. 

From an arbitrary starting point, a random step is taken and the cost function is evalu- 
ated. If the step has improved the cost function (a success), the present parameter val- 
ues are used as a new starting point and a new step is taken. If the cost function was not 
improved (a failure), the original parameter values are retained, and a new step is taken. 
Thus a "walk” over the contour is generated. 
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For explanation of the basic algorithm in more detail, attention is directed to fig- 
ure 2(b). As illustrated in the figure, there are nine possible directions that can be taken 
from a given starting point in the two-parameter plane. These nine directions are 
obtained by plus, minus, and zero perturbations of each parameter. (Note that one direc- 
tion is no movement at all.) The eight directions that involve movement are stored in an 
array, which the algorithm will sample randomly without replacement. The array is 
sampled until a success is obtained or all directions are exhausted, at which point the 
array is reinitialized. If failures should occur in all possible directions from a given 
point, a reduction in step size would be necessary in order to continue iteration. The 
algorithm is stopped once the step size becomes smaller than the desired detectable dif- 
ference (see eq. (13)) between different values of the cost function. 

In conjunction with the basic algorithm, two features are available to speed conver- 
gence. The first feature is a strategy implemented at the option of the designer by the 
algorithm and simply requires that after each success the next try be made in the same 
direction. The second feature involves interaction of the beam designer with the algo- 
rithm through the CRT displays and will be described later. 

Parameter Constraints 

The parameter constraints of equations (8) are implemented by simply not allowing 
the algorithm to perturb parameters across boundaries. Thus these constraints affect 
the problem only in that the parameters must always be within their feasible ranges. The 
focal-point constraints, equations (9), are implemented by perturbing parameters, calcu- 
lating the focal points, and sensing for violations. Should a violation occur, the perturba- 
tion direction is deemed a failure, the equations of motion are not evaluated, and a new 
direction is chosen by the algorithm. It should be obvious that the addition of the focal- 
point constraints greatly speeds convergence, since the evaluation of the equations of 
motion is bypassed for all parameter sets which do not satisfy the constraints. 


Statistical Test for Success 


Assuming equal variances, a pooled t test for equality of means, as described in 
reference 11, is used to determine the success or failure of a perturbation direction. The 
null hypothesis was chosen to assure control over the probability of accepting a failure as 
a success (a 5-percent chance of occurrence was selected). The test is shown as follows: 
Hq: Cp S c (Failure) H^: Cp > C (Success) 


*cal “ 


Cp - C 


l (JL + l\ 

( n p - ijop + (n - 1)0-2 

' \ n P n / 

np + n - 2 


(ID 


Reject Hq if t ca i > t. 05 (118) = 1.66 
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where rip + n - 2 = 118 degrees of freedom. Here, Cp is the value of the mean cost 
at the last success; n p is the sample size (60 rays/beam) at the last success ; d 2 is 

the variance at the last success; and C, n, and &2 are the mean cost, sample size, 
and variance, respectively, for the present parameter values being tested. It should be 
noted that an F test was used periodically to validate the assumption of equality of 
variances, as shown in the following equations: 


H 0 : o^ = cr2 


H a : / a2 


r 


F . = — where < 
cal b N 


A = ma x(o^,ct2 


B = min^5p,o^j 
Reject H 0 if F cal > F >95 (60,60) = 1.67 


(12) 


The null hypothesis was accepted in all cases, except for a few in which the step size of 
the perturbations was large. 

The sample size (60 rays/beam) was determined by using the iterative method 
described in reference 12, pages 154-155, 


> 2 ( t o + h) 2 ^ 2 
6 2 

where 


(13) 


t Q t value associated with type I error; t 0 = t qs(1 18) = 1.66 

tj t value associated with type n error (95-percent power for the chosen detectable 

difference); tj = t i 0 (118) = 1.29 

$2 estimate of variance; d 2 “0.01 

6 detectable difference desired; 5 = 0.055 

GRAPHIC DISPLAYS AND INTERACTIVE CAPABILITIES 

Extensive use of the CRT displays is an integral part of program operation. In 
addition to the algorithm displays, the following problem -related displays are available 
to the analyst: 
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(1) Phase plots for both the x and z variables at the entrance to the transport 

system and at the exit of each quadrupole 

(2) An xz-plane plot at the same locations 

(3) Beam traces, both horizontal and vertical views, including the location and size 

of each quadrupole, slit placement in the system, and the horizontal and verti- 
cal focal points 

Thus the analyst has complete information of the transport system at important points 
along the system. If the designer desires, based on his knowledge of transport systems, 
algorithm operation may be bypassed, and with use of the displays, analog matching may 
be carried out. This "hands-on" capability allows the designer to look at areas in the 
parameter contour which he may deem interesting. 

The algorithm displays are available to possibly speed convergence to an extremum 
of the problem in the following manner: The designer may, during algorithm operation, 
display on the CRT a history of the cost-function behavior versus the number of tries 
(both successes and failures) and also histories of the value of each parameter versus the 
number of successes. Based on this information, the designer may be able to determine 
the sensitivities of certain parameters in the current region of operation of the cost 
domain. If the parameter appears sensitive, the gain, or perturbation size of that param- 
eter, may be increased. On the other hand, if a certain value is required of a parameter 
in order to obtain a success or if the parameter appears insensitive (varies randomly 
about some mean), it may be made a constant. In this situation, the number of possible 
directions is greatly reduced, as evidenced by table 1, and the array of directions is 
regenerated for the reduced number of parameters. Therefore, increasing the gain of 
sensitive parameters and eliminating insensitive parameters should then speed conver- 
gence to an optimal set of parameters. 


RESULTS 

The major fields of interest in beam-transport design are transporting and trans- 
forming either a parallel source or a point source, at the exit of the accelerator, into a 
parallel distribution or a point distribution at a desired location. Thus, four cases arise: 
parallel to parallel, parallel to point, point to parallel, and point to point. All four are 
presented to illustrate the versatility of the method. 

Figure 3 shows a typical sampling from the initial distribution used in all parallel- 
source examples except the parallel-to-point case. Although 60 points were used in each 
algorithm iteration, only 30 are shown in order to provide clarity in the plotting. The 
distribution used to model the beam-entrance conditions was fixed along the envelope of 
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the phase-plane area, since interior points in the phase -plane area remain in the interior 
throughout the transport system. 

Four studies were made under the parallel-to-parallel case, all of which used the 
equations 
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as the cost function defining the desired end conditions. (The desired end conditions are 
represented as an ellipse in the appropriate phase-plane figures.) The four studies are 
(1) best analog results, (2) best digital results, (3) addition of chromatic aberration to the 
digital problem, and (4) addition of slits to the digital problem. 

Table 2 together with figures 4 to 7 shows a comparison of typical results at the 
end of the last quadrupole obtained from the analog program and the digital program using 
the analog-obtained parameter values as a starting point for the digital program. As 
exemplified in table 2, the digital program in all cases compared made changes in the 
parameters which lowered the value of the cost function. Comparison of figures 5 and 7 
reveals the visual-resolution problem mentioned previously. Figures 8 and 9 are results 
obtained by adding 5 percent chromatic aberration to the digital system of figures 6 and 7. 
Table 2 also lists pertinent information for this case. During this study it was noted that 
in most cases the introduction of chromatic aberration did not affect the convergence of 
the algorithm. 

Table 3 and figures 10 to 13 illustrate an example of the introduction of a slit in one 
plane to achieve some desired end conditions which cannot be met by the system without 
slits. In the example the requirement to eliminate all slopes greater than a certain value 
in the vertical plane, as shown by the points outside the ellipse in the lower portion of fig- 
ure 10, has been placed on the system. Assuming it has been determined that the system, 
subject to parameter constraints, cannot meet the requirement, the problem is solved in 


the following manner: First, the additional constraint that the focal point of the first 
doublet in the vertical plane lie within a certain area between the two doublets is intro- 
duced into the system (eq. (9b)). The algorithm is then used to determine a set of system 
parameters satisfying all constraints which yield end conditions best fitting the imposed 
requirement. At this time a slit is placed at the vertical focal point of the first doublet, 
and the slit opening is varied until the imposed requirement is met and the desired end 
conditions are achieved (fig. 12). 

The last parallel -to -parallel case is a comparison of computer-generated results 
with those obtained from an actual experiment conducted at the NASA Space Radiation 
Effects Laboratory. Figure 14(a) reflects the computer -generated results while fig- 
ure 14(b) depicts the results obtained from the actual experiment. (Fig. 14(b) shows the 
proton beam passing through a spark chamber after it has exited the beam transport 
system.) 

In the point-to -parallel study, equations (14) are still used as the cost function since 
the desired end conditions are the same as in the parallel -to -parallel case. Therefore, 
the only change required to enable the program to go from the parallel-to -parallel case 
to the point-to -parallel case is a change in initial conditions. The algorithm will then 
determine a new set of parameters which will best meet desired end conditions. Fig- 
ure 15 shows a typical sampling from the initial distribution used to simulate a point 
source. Figure 16 reflects the end conditions obtained, and table 4 lists the set of system 
parameters obtained by the program in addition to values of other variables for the point- 
to -parallel case. Figure 17 gives a good visual illustration of the point-to -parallel case. 

The parallel -to -point and point-to-point cases are obtained by using the appropriate 
initial conditions together with a cost function defining the desired end conditions as a 
point. The following equations represent the cost function used in these studies: 
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Figure 18 is a typical sampling of the initial distribution used for the parallel-to -point 
case. Figure 19 and table 5 show the results obtained from the program. It should be 
noted that the distance D4 between the end point and the beginning of the first magnet 
is an additional parameter to be determined in both the parallel-to -point and point-to- 
point cases. 

The final case, point to point, is shown in figure 20 with its pertinent data given in 
table 6. The initial conditions used in this case are the same as for the point-to -parallel 
case shown in figure 15. 


CONCLUDING REMARKS 

A modified creeping-random-search algorithm has been applied successfully to 
solving a particular beam -transport design problem. The method of application that 
evolved in the process of solving this problem possesses the advantages of algorithm 
operation with analyst interaction or if desired, analog matching with digital accuracy, 
and feasible solutions in the presence of both position- variable and parameter constraints. 
In addition, realistic modeling of chromatic aberration has been included in the mathemat- 
ical model as affecting individual members of the beam, rather than the beam as an entity. 

The ability of the method to handle all major cases has been demonstrated using a 
particular mathematical model, although flexibility in the model is possible. Results 
presented have shown that this method not only yields better solutions than those of a 
standard analog program, but also yields solutions for systems, subject to constraints, 
for which prior transport optimization programs give infeasible solutions. Also illus- 
trated by the results are the introduction of chromatic aberration and the use of a slit to 
obtain desired end conditions. The results demonstrate that good agreement was obtained 
between the results predicted by the program and actual experimental results obtained 
with the 600 MeV cyclotron at the NASA Space Radiation Effects Laboratory. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., November 12, 1971. 
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TABLE 1.- VARIATION OF NUMBER OF POSSIBLE DIRECTIONS 
WITH NUMBER OF PARAMETERS 


Number of 
parameters 

Possible directions 

w 

i 

2 

2 

8 

3 

26 

4 

80 

5 

242 

6 

728 

7 

2186 


3 n - 1 possible directions. 
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TABLE 2.- PARALLEL-TO-PARALLEL CASE 


Parameter 

Analog 

Digital 

No aberration 

5% aberration 

Kj, m -1 

0.9110 

1.1525 

1.1525 

K 2 , m -1 

-0.9398 

-1.5665 

-1.5665 

K 3 , m -1 

0.8812 

0.9125 

0.9125 

K 4 , m-l 

-0.6596 

-0.8250 

-0.8250 

^ 2 , m 

3.21 

4.0125 

4.0125 

C 

3.2747 

2.2222 

3.2117 

F ( s f) 

1.4965 X 10-3 

1.0562 x 10-3 

1.2758 X 10-3 

z '( s f) 

1.59914 x 10"3 

V 

1.0949 x 10"3 

1.1208 x 10-3 

J 

h> m 




0.7176 


m 


0.7176 


h> m 


0.6858 


Z 4 , m 


0.6858 


Dl> m 


0.1969 


E>3> m 


0.2286 


N, percent 


0 


w dif 


1.0 


W N 


0 
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TABLE 3.- SLIT IN VERTICAL PLANE 


Parallel-to -parallel case; no chromatic aberration 


Parameter 

Without slit 

With slit 

3 

1 

1.165 

1.165 

K 2 , m-1 

-1.554 

-1.554 

K 3 , m -1 

0.938 

0.938 

K 4 , m“l 

-0.837 

-0.837 

D2> m 

3.747 

3.747 

C 

2.467 

1.967 

F ( s f) 

1.119 X IO - 3 

1.027 x 10 -3 

z '( Sf ) 

1.276 X IO - 3 

6.566 X IO - 4 

N, percent 

0 

63.33 


x 

J 




* 1 » m 

0.7176 

^ 2 j m 

0.7176 

^3> m 

0.6858 

Z 4 , m 

0.6858 

Di, m 

0.1969 

D3> m 

0.2286 

Wdif 

1.0 

r,m 

0.3 

x sl> m 

0.045 
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TABLE 4.- POINT-TO-PARALLEL CASE 


[No chromatic aberration; no focal-point constraints] 


Parameter 

Value 

K l> m -1 

1.265625 

K 2 , m -1 

-1.065625 

K 3 , m -1 

1.118750 

K 4 , m~l 

-0.934375 

D 2 > m 

4.28750 

F ( s f) 

1.1102 X 10-4 

^(Sf) 

2.2458 x lO - 4 

C 

3.14786 

^ 1 , m 

0.8128 

ll, m 

0.8128 

h> m 

0.5588 

m 

0.5588 

Dl> m 

0.1016 

D3, m 

0.0508 

Wdif 

1.0 

W N 

0 

x(0), m 

0.001 

x'( 0 ) 

0.01 

z( 0 ), m 

0.001 

z'( 0 ) 

0.01 


23 


I 


TABLE 5.- PARALLEL-TO- POINT CASE 


[No chromatic aberration; no focal-point constraints] 



Value 
-0.2340 
0.4180 
-1.0650 
1.4468 
4.7875 
9.0063 
4.2084 x 10-7 
0.5588 
0.5588 
0.8128 
0.8128 
0.0508 
0.1016 
1.0 
0 

0.055 

0.000235 

0.055 


0.000235 



TABLE 6.- POINT-TO-POINT CASE 


| No chromatic aberration; no focal-point constraints] 


Parameter 

Value 

V 

h* 

3 

1 

H* 

0.4000 

K 2 , m -1 

1.3000 

K 3 , m ' 1 

- 1.2000 

i-H 

1 

a 

1.2499 

E> 2 , m 

4.7000 

D 4 , m 

11.0656 

C 

6.453 X 10 ‘ 6 

h> m 

0.7176 

h, m 

0.7176 

lz, m 

0.6858 

£4, m 

0.6858 

Dl» m 

0.1969 

D 3 , m 

0.2286 

Wdif 

1.0 

W N 

0 

x( 0 ), m 

0.001 

x'( 0 ) 

0.01 

z( 0 ), m 

0.001 

z'( 0 ) 

0.01 


25 



Quadrupole magnetic fields 



* is a variable in the parallel-to-point and point-to-point cases 



Quadrupole axis system 

Figure 1.- Beam transport system. 
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Figure 4.- End conditions for analog case. Kj = 0.9110 m~l; 
K 2 = -0.9398 m-1; K 3 = 0.8812 m"l; K 4 = -0.6596 m~l; 
and Do = 3.21 m. 











z 


Figure 8 .- End conditions for digital case with 5% chromatic aberration. 
K x = 1.1525 m-1; K 2 = -1.5665 m-1; K 3 = 0.9125 m-1; 

K 4 = -0.8250 m _ l; and D 2 = 4.0125 m. 
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Figure 9.- Trace of digital case with 5% chromatic aberration. Kj = 1.1525 m - l 
K 2 = -1.5665 m-1; Kg = 0.9125 m-l; K 4 = -0.8250 m" 1 ; and D 2 = 4.0125 m 


















Figure 15.- Distribution for point-to -parallel case initial conditions 





BOTTOM VIEW 



SIDE VIEW 

Figure 17.- Trace of point-to -parallel case. Kj = 1.265625 m - l; 
K 2 = -1.065625 m-1; K 3 = 1.118750 m-1; K 4 = -0.934375 m-1; 
and D 2 = 4.28750 m. 
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SIDE VIEW 


Figure 19.- Trace of parallel*- to -point case. 
K 2 = 0.4180 m-1; K 3 = -1.0650 m“l; K 4 
D 2 = 4.7875 m; and D 4 = 9.0063 m. 


Kj - -0.2340 m-1; 
1.4468 m“l; 








m 







m 




BOTTOh VIEW 



SIDE VIEW 

Figure 20.- Trace of point-to-point case. K 4 = 0.4000 m~l; 
K 2 = 1.3000 m-1; K 3 = -1.2000 m-1; K 4 = 1.2499 m-1; 
D 2 = 4.7000 m; and D 4 = 11.0656 m. 
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